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Abstract 

We develop analytical tools and numerical methods for time evolving the total density matrix of the 
finite-size Anderson model. The model is composed of two finite metal grains, each prepared in canonical 
states of differing chemical potential and connected through a single electronic level (quantum dot or im- 
purity). Coulomb interactions are either excluded all together, or allowed on the dot only. We extend this 
basic model to emulate decohcrring and inelastic scattering processes for the dot electrons with the probe 
technique. Three methods, originally developed to treat impurity dynamics, are augmented to yield global 
system dynamics: the quantum Langevin equation method, the well known fcrmionic trace formula, and an 
iterative path integral approach. The latter accommodates interactions on the dot in a numerically exact 
fashion. We apply the developed techniques to two open topics in nonequilibrium many-body physics: (i) We 
explore the role of many-body electron-electron repulsion effects on the dynamics of the system. Results, 
obtained using exact path integral simulations, are compared to mean-field quantum Langevin equation 
predictions, (ii) We analyze aspects of quantum equilibration and thcrmalization in large quantum systems 
using the probe technique, mimicking clastic-dephasing effects and inelastic interactions on the dot. Here, 
unitary simulations based on the fermionic trace formula are accompanied by quantum Langevin equation 
calculations. 
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I. INTRODUCTION 



There has been recently a great deal of interest in simulatingthe real-time dynamics of quantum 
systems, open or closed, prepared in a nonequilibrium state These investigations have been 
spurred by recent experimental breakthroughs in the ability to watch out-of-equilibrium dynamics, 
for example, in cold atomic gases [2], or in on-chip superconducting circuits [jj]. This endeavor 
is fundamentally important for resolving basic issues in quantum dynamics, and in particular, 
for understanding equilibration and thermalization in quantum systems The nonequilibrium 
dynamics of the eminent Anderson model J], composed of a single electronic level (quantum dot) 
coupled to two metals, has in particular been of great interest. This is because it is perhaps the 
simplest platform for probing both equilibrium and out-of-equilibrium physics in a many-body 
system. The model is integrable, even when electron-electron (e-e) repulsion effects are accounted 
for on the dot, and its integrability has been exploited for resolving its transport behavior 

A central strategy in most analytic and numerical tools devoted to the Anderson model, and 
impurity systems at large, is the separation of the total system into a subsystem (dot) and the 
environment (metals, referred to as reservoirs). The latter are typically assumed to be infinite- 
dissipative and are maintained in one of the canonical ensembles of statistical mechanics. This 
assumption allows one to treat the effect of the reservoirs on the subsystem within a self-energy 
term. However, once the reservoirs are traced out, one cannot describe their explicit dynamics. 
Among the numerical approaches developed along these lines we list the time-dependent numerical 
renormalization-group method [fj|, real-time diagrammatic Monte Carlo techniques [?J and path 



integral approaches [8|, M\. These methods place focus on quantities such as the dot occupancy, 
transmission probability, conductance, current, noise, and correlations on the impurity. The dy- 
namics of the total system, including the electron reservoirs, has not yet been explored since general 
tools for simulating the overall dynamics in a system-bath scenario are still missing. 

The current work develops analytical and numerical treatments of global system evolution based 
on established impurity dynamics techniques. These tools allow investigation of the roles of e-e 
interactions and decoherence and dissipation effects on nonequilibrium reservoirs dynamics. We 
focus on the finite-size Anderson model composed of two metallic grains weakly coupled through 
a single electronic level. We refer to the metal grains, each composed of N ~ 100 — 500 electronic 
states and n ~ 50 — 200 electrons as "reservoirs" alluding to their high density of states (DOS). 
This large DOS allows the reservoirs' effect on the dot (subsystem) to be absorbed into a p ositive 



real self-energy function lending to a quantum Langevin equation (QLE) description 



10h12|] , as we 
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FIG. 1: Two metallic grains (reservoirs) separately prepared in a grand canonical- diagonal state. At to the 
reservoirs are put into contact through a single electronic state. We study three variants of this systems: a 
case without interactions, (a) allowing for electron repulsion effects on the dot only, and (b) accommodating 
decoherring and inelastic effects on the dot, by coupling its electrons to a G reservoir, serving as a dephasing 
or a voltage probe. 

explain below. A schematic representation is presented in Fig. [TJ We are interested in following the 
real-time dynamics of both the dot and the reservoirs degrees of freedom. As an initial condition, we 
assume that each reservoir is prepared in a distinct Gibbs-like grand canonical state at a different 
chemical potential but at the same temperature. 

In the absence of dephasing and inelastic effects, the dynamics of the total density matrix 
is followed by extending three approaches: (i) the quantum Langevin equation method 1CH12||. 



adopted here both in the noninteracting limit and in the mean-field (MF) regime, (ii) fermionic 



trace formula [13J], used here for simulating the exact dynamics of the noninteracting model, and 

nn 

(iii) an influence functional path integral method [141 . |15| ] , employed to treat interactions beyond the 
perturbative regime. In the latter half of the paper, these techniques are used to study reservoir 
population evolution both without and with Coulomb repulsion effects on the dot, and in the 
presence of emulated dephasing and inelastic scattering effects. 

While Coulomb interactions are explicitly introduced here, the inclusion of dephasing and in- 
elastic effects warrants further discussion. The origin of such processes are many-body interactions 
in the system, e.g., electron phonon coupling. Since an explicit and exact inclusion of these in- 



teractions is extremely challenging 
their stand 



161-4191] . phenomenological techniques have been developed in 
2CH22l]. In the case of elastic decoherring processes the technique is referred to as 
a "dephasing probe". In the case of inelastic scattering processes, it is referred to as a "voltage 
probe". These probes are electron reservoirs, prepared such that, there is no either energy resolved 



or total net electron flow from the L-dot-R system towards these probes. For a scheme of this 
model, see Fig. QJb). It should be noted that elastic-decoherring processes or inelastic effects are 
only emulated here by t he p robes. The overall dynamics can be still simulated using the unitary 
trace formula technique [23J. We also extend the QLE method to include a probe, and time evolve 



the system. Since our calculations provide the real-time dynamics of the full density matrix (DM 
the process of equilibration and thermalization in a finite quantum system can now be studied 
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261 ] . Particularly, we find that when only decoherence effects are allowed, the system approaches a 
non-canonical equilibrium state. In contrast, when inelastic processes are included, the reservoirs 
relax towards a common Gibbs-like state. 

The paper is organized as follows. In Sec. [II] we present the finite-size closed Anderson model 
and outline the implementation of the probes. In Section [TTT1 we present our developed numerical 
and analytic treatments of the density matrix dynamics: First, we extend the standard quantum 
Langevin equation approach to include reservoirs dynamics. The method can treat both the nonin- 
teracting model A the case with interactions, only at the level of mean-field (Hartree) theory 
(|III A 2j) . and the probe model (jlll A 31) . Second, we present the fermionic trace formula, useful for 
studying the Anderson model without interactions and with implemented dephasing and inelastic 
effects in Sec. IIIIBI The third method, presented in Sec. IIII C\ is an influence-functional path 
integral approach 1J) 3j- This non-perturbative tool can treat the model with interactions in a 
numerically exact manner. Applications are included in Sec. IIV1 The effects of Coulomb repulsion 
effects on the dot are studied using mean-field QLE and the path integral technique in Sec. IIV Al 
In Sec. IIV Bl quantum equilibration and thermalization is investigated using the probe technique. 
In this case, the total density matrix is resolved using the QLE method and the fermionic trace 
formula. The paper is summarized, along with an outlook, in Sec. [Vj 



II. MODEL 



The closed-system Anderson model consists two metal grains, v = L, R, including (each) a 
collection of N u dense electronic levels initially populated by noninteracting electrons up to the 
chemical potential fj, u , at temperature T = f3~ 1 . The two baths couple only through their (weak) 
hybridization with a single level quantum dot. Work presented in this study concerns three variants 
of the model. The simplest version is the "noninteracting case" , where electron-electron repulsion 
effects and any decoherring and relaxation mechanisms are excluded. The second case, the "inter- 
acting model", allows for e-e interactions on the dot only. The third model variant, the "probe 
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model", phenomenologically contains elastic decoherring and inelastic scattering processes on the 
dot, using the probe technique and excluding e-e repulsion effects. This model is discussed in detail 
in Sec. IIVB1 where it is applied in the context of quantum equilibration. 

A. The interacting model 

In the absence of decoherence and dissipation effects the interacting Hamiltonian takes the form 

H = H L + H R + H W + V L + V R , (1) 

where Hl,r,w represents the Hamiltonian for the left reservoir, right reservoir and the dot, respec- 
tively. The term V u denotes the coupling of the dot to the v reservoir, 

Hl = y"]^4(T ( V> H R = yZ e r c t,a c r,a 
I, a r,a 

V L = ^2 V l C d,o C l,° + h - C - V R = ^2 V r C d,a C r,« + h.C. 
la ra 

Hw = X! tdc\ a c d)(T + Un d ^n d ^. (2) 
u 

Here, C)- a (k = l,r,d) are fermionic operators of the left reservoir, I € L, right reservoir, r € R 
and the dot (d). The symbol a stands for the spin state (f or J,) and U accounts for the onsite 
repulsion energy. We assume that v\ and v r are real numbers and that the Hamiltonians of the leads 
are diagonal in momentum basis and define the hybridization r„(e) = ir Ylkev v k^( e ~ e k)i taken 
in practice to be energy independent. The Hamiltonian ([2} disregards magnetic fields, yielding 
spin-degenerate energy levels, thus it is sufficient to consider observables for one spin species. We 
note that the noninteracting case arises simply from the suppression of U. 

Our objective in this paper is to calculate the time evolution of the expectation values of all 
two-body operators in the system (k,j = l,r,d) 

Pk,j(t) = <4(t)c,-(t)> = Tr[p(t )cl(t) Cj (t)}, (3) 

written here in the Heisenberg representation with p(io) = Pd <8> Pl® Pr representing the factorized 
time-zero density matrix of the system, and with the trace performed over all degrees of freedom. 
We suppress the spin degree of freedom in the density matrix since its elements are identical for the 
two spin configurations. As an initial condition, we take the dot to be empty and the reservoirs' 
DM to be diagonal, 

(C d (t )c d (t )} = 0, (4(t )Q(t )) = = fl, (4(t )Cr(t )) = /fl(Cr) = fr, (4) 
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with the population /r,,it(e) = [e^ e vl,r) _|_ i] As a convention, we use the symmetric chemical 
potential bias [Il = > 0. 

B. The probe model 

The Anderson probe model, a variant of the basic model, Eq. ([2]), can emulate memory loss and 
energy redistribution in a quantum system without explicitly introducing many-body interactions 



2CB22I]. The probe technique has been of extensive use in meso scop ic physics, for describing the 
disappearance of quantum effects in transport [21[, dissipation 20(], and equilibration dynamics 
[27! . Recent advances include a full-counting statistics analysis of the probe model 



and an 



29j. We model here either a dephasing probe, 



extension of the probe technique to the AC regime 
allowing for quasi-elastic decoherence processes, or a voltage probe, where inelastic effects are 
further mimicked. In both cases we suppress electron-electron interaction effects in the system. 

As we explain next, in our study the "probe" terminology refers to a setup slightly different 
from the conventional one. The standard construction refers to an open system scenario, where 
the probe practically performs which-path experiments through repetitive measurements of the 



system [30f|. In contrast, in our picture the probe is a finite-closed quantum system, only initialized 
with a certain-special distribution. After its preparation, the probe, similarly to other parts of the 
system, is left undisturbed. Thus, we can use exact unitary approaches and simulate the dynamics 
of the total system. While this picture abuses to some extent the standard notion of a "probe" , 
we maintain this terminology here since practically our implemented probe acts like a proper one, 
inducing phase loss or/and energy reorganization in the system. 

We introduce a probe into the model by adding an additional Fermi-sea reservoir, denoted by 
the letter G, to the Hamiltonian (J2]), again discarding the spin degree of freedom, 

H P = H + H G + V G , (5) 

where 

H G = ^ e g c l c 9> V G = ^2 v 9 c l Cd + h ' c - ( 6 ) 
9 9 

We naturally define the hybridization T^e) = ^Sg^g^C 6 — e a)' anc ^ take it as a constant. As 
always, our objective here is the resolution of all system expectation values of two-body operators 
(k,j = l,r,d), pkj(t). As initial conditions we assume Eq. (j4|), where the G bath initial condition 
is set according to the particular probe condition, explained below. 
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Voltage probe. Inelastic scattering effects of electrons on the dot are effectively included by 
implementation of a voltage probe. The probe has a canonical distribution, /g(e) = [e^ e-/iG )+l]~ 1 , 
and its chemical potential [ic is set such that the net charge current from the dot to the G unit 
vanishes for all times 

i G=J t Y,( C l C 9)=°- ( 7 ) 
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With the motivation to explore situations beyond the linear response regime 
numerically, by employing the Newton-Raphson method 321 ] . 



311 ]. we retrieve \xg 



(m+l) _ (m) . (.im)\i-i / (m)s ( q\ 

IIq is the initial guess, i G denotes the first derivative with respect to fic- I n principle, one should 
adjust fj>G throughout the simulation, to eliminate population leakage from the L-dot-i? system 
into G. However, we have found in our simulations that the G bath has lawfully behaved as a 
probe once we determined \iq from the steady-state limit using the following analytic expression 
for the charge current 

. f , 2r G T R [f G (e) - f R (e)] + T L [f G {e) - f L {e)} 

ic = J «c(e)de, (9) 

with r = + T/j + Tq 33J. The lower and upper integration limits are determined by the band 
simulated. Substituting Eq. ([9]) in Eq. (J7|), a voltage probe condition is set by demanding /g(c) 
to fulfill the relation 

(e-( d p+P T L +T R J (e- ei ) 2 +r 2 1 ' 

Dephasing probe. Implementation of the dephasing probe, fabricating elastic decoherence, necessi- 
tates the stronger requirement «g( £ ) = 0? i- e -> t ne charge current at a given energy should vanish. 
Using the steady-state behavior ([9]), we obtain a non- Fermi distribution 

f / a rg/g(e) + T L f L (e) 

fc(e) = • (11) 

We emphasize that \xq or fc have been determined here in the steady-state limit, assuming fixed 
chemical potentials for the L and R baths. Indeed, at short time, T^^t < 2, before a (quasi) 
steady-state sets in, we find that iq ^ 0. However, we have confirmed numerically that beyond 
this time throughout all our simulations \ic/iL,R\ < 10 -4 , thus the G reservoir plays the role of a 
proper probe. 

Three different approaches for the calculation of the full DM are described in Sections IIII Al 
iHJBl and lHTCl Applications are included in Sec. IIVI 



III. METHODS 



A. Quantum Langevin Equation 

The dynamics of the Anderson model in the absence of interactions (U = 0), with interactions 
at the mean-field level, or with a probe, is described here within a quantum Langevin equation 
framework [l(3|. The basis of our method has been used in the past to follow the dot evolution 
or the charge and energy currents in the system 
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12 1 . Here, we show results for the full 



DM. We begin our analysis with the trivial treatment of the impurity (dot) and review the steps 
involved. This review helps highlight underling approximations and establishes limits for the 
method's applicability. 

1. Noninteracting case (U = 0) 

In the Heisenberg representation the fermionic operators satisfy the following equations of mo- 
tion (EOM), 



T 

r 



Q = -ie L ci - iv t c d 

c r = —ie r c r — iv r Cd (12) 
Formal integration of the reservoirs EOM yields, e.g. at the L end, 

j(t) = e-*i(*-*0q(i )-iw, [ dre- ie ^c d (T)dT. (13) 



Cl{ 



We substitute Eq. (fT3|) . and the analogous expression for c r (t), into the dot EOM [Eq. (fT2|) ]. and 
retrieve 



c d = -^Q-z^^e-^^^cKto)-^^-"^*-* ^,.^) 

I r 

~ f dT^v?e- lei(t - T) Cd(T)- f drY.vle-^-^Cdir). (14) 

J to i J to r 

In this (exact) equation the second and third terms are interpreted as "noise" [lo| . 

r? L (i) = J^e-^-^ciCto) 
I 

r! R (t) ee J2 v re~ ierit ~ to) Cr(to). (15) 
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The last two terms in Eq. (|14p can be reduced, each, into decay terms, further inducing an 
energy shift of the dot energy, absorbed into the definition of e^. This is justified by following two 
assumptions: (i) The hybridization r„(e) = tt *}2 kei/ v k 5(e — ej.) may be taken as a positive real 
self-energy function [hJ, and (ii) the dot dynamics is slow relative to the reservoirs' evolution. We 
now explain the steps involved. First, we define a new operator for the dot, by absorbing its fast 
oscillatory behavior, Cd(t) = Cde tedt . Its EOM is 

't rt 



Cd 



-i[r} L (t) + 7] R (t)]e 



v 2 e -ie ld (t-r) 5d ^ - f drJ2 
J to r 



y 2 -ie rd (t-r) 



Cd(r) (16) 



where e k j = £fc — £j. We then change variables, x = t — r, and make the assumption that the dot 
evolution (now missing the fast phase oscillation) is slow with respect to other time scales in the 
system, Cd(t — x) ~ Cd(t). This results in 



2 e -ieid(t-r)~ 



'to l 

If the time is long, t S> to, integration gives 



t-to 



dx 



e - itldX dx = T L (e d ) - 2ivf lim 



t— >oo 



D L (e) 



sin 2 (e 



e d )t 



de 



with Di{e) = Y2i <5( e ~~ e i) as the density of states of the L metal, taken as flat here. We also take 
the interaction parameters vi to be independent of the I index. The imaginary term introduces an 
energy shift, which can be absorbed into the definition of e^. It diminishes when the density of 
states does not depend on energy (the case used later), and when the bandwidth is large enough. 
In our numerical calculations we have used a finite bandwidth with a cutoff D = ±1, introducing 
a small correction to e^. We return to Eq. (JT3J) and conclude that it obeys the quantum Langevin 
equation 



Cd 



-ie d c d - irj, L {t) - ir) R {t) - T(e d )c d (t), 



(17) 



with T(e) = r^(e) + T^(e). The dynamics of the dot occupation, (n^(t)), can be reached by a 
formal integration of Eq. (|17p , 



c d (t) = c d (t )e ( -^- r)(i - io) 
to provide the standard expression [341 ] 



lie 

h 



(-ie d -T){t-T) { L 



[rf(T) + r,"(T)]dT, 



(18) 



(n d (t)) = (c d (t)c d (t)) = Y, 



\vk\ 2 fk 
- e 2 +r 2 

k=l,r dk 



1 + e 



-2T(t-to) 



2e 



-r(t-to) 



cos[€ d k(t - t )] . (19) 
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This derivation relied on the initial conditions (j3]). The summation runs over all the reservoirs 



degrees of freedom [35( . We now use Eq. (fT3"j) and its analogous expression for c r (t) , together with 
Eq. (|18p . and derive analytical expressions for other quadratic expectation values, {c\(t)cj(t)}, 
k,j = l,r,d. These results are valid as long as one can faithfully rely on Eq. (|17p . In what follows 
we take to = to simplify our notation. The reservoir-dot coherence can be obtained analytically, 



p dtl (t) = (cl(t) Cl (t)) = B X + B 2 . (20) 
Here, B\ includes contributions from the L side only, 



D ivifi 

-D1 



1 - e 



-t(F-ie dl ) 



(21) 



T - ie d i 

B2 includes electron transmission pathways from the L side, through the dot, to the K = L,R 
grain, 



b 2 = -wi 



v lfk 



r 2 + e 



2 



{I _ pitklt e~ 2Tt — e _i T-* e di) g-^r-iefed) _ p-i^lkt e -t{T-ie dk ) _ e -t(r+ie M ) ~) 
— + : 7; : p : >■ (22) 

Using Eq. (I20p . we derive an expression for the charge current at the L contact, 

= |E^WQW) = -2S^^(4(t)Q(t) 

2 rLrj; \ - fi - /g _ 2r L _ rt \ - 1 

x J e" r * (/,r L + / r r fl ) - {2hT L + 2/ r r fl - r/,) cos (e w t) - /,e w sin (e w *) I . (23) 



Here S stands for the imaginary part. An analogous expression can be written for in(t). Ref. 
includes a Green's function based derivation for the time dependent current in the symmetric limit 
(r^ = r^j). This derivation results in a surplus nonphysical term at the initial time. We now turn 
our attention to the reservoirs' states population. Using Eq. (fl3j) . we find that it is given by three 
contributions, 

p(e,) = (4(i)q(*)) = (cf(to)q(io)) 

+ ivie- itl{t - to) [ e iei{t ~ T) {c ] d {T)ci{t G ))dT + c.c. 
Jt 

+ ^ 2 /" rrfr 1 dr 2 (c^(r 1 ) Cd (r 2 ))e^^)e-^( t - T " 2 ). (24) 
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The two-times correlation functions can be obtained from Eqs. (|13p and (|18j) without additional 
approximations, and the explicit expressions are given in Appendix A. Similarly, closed analytic 
expressions can be written for inter and intra-reservoir coherences, e.g. pik(t), k 6 L,R, see 
Appendix A. Eqs. ([19]), ([20]) . (|24p. (|A3|) and (|A4|) . and the analogous /2-bath expressions, form 
the time-dependent full density matrix of the system. 

Timescales. We now comment on the applicability of the QLE approach. Given infinite reser- 
voirs, a current carrying steady-state behavior develops, and the dot occupation, as well as the 
charge current, reach a fixed value after a short time, t% > 2/T (see for example Fig. [8]). However, 
since the reservoirs are finite in the present treatment, recurrence effects should eventually man- 
ifest in our system. These effects cannot be handled by the QLE technique since an irreversible 
behavior has been assumed for the dot, as Eq. (|17p breaks unitary evolution. The technique can 
still excellently reproduce the exact dynamics in the so called "quasi steady-state" (QSS) region, 
up to ~ 2ttN/D 36|]. Here iV is the number of electronic states in each bath and D is the 
band cutoff. Within this time, the dot occupation and the charge current are constant, similar to 
a real steady-state situation. Around the time the dot occupation should begin to vary, showing 
(partial) recurrence behavior, and the QSS limit breaks down. 

Fig. [2] clarifies this timescale issue. The left panel displays the dot occupation as a function 
of time using either the QLE approach (full line) or an exact method (dashed line), described 
in Sec. IIIIBI Results agree up to t ~ ~ 630, and deviations before this time are due to the 
finite band used in QLE, while neglecting the energy correction in Eq. (|17p . Around the time 
exact simulations show a partial reversal of the dot occupation, while QLE still produces the QSS 
value. At a later time, t ^> r<2, the QLE data diverges. Panel (b) presents the bath occupation for 
two selected energies, and we find that nonphysical values, such as a population exceeding unity, 
can be obtained with QLE when t > t&. Thus, the QLE method can be used within the interval 
t < only, to be consistent with its underlying assumptions. However, interesting nonequilibrium 
physics takes place within this window, thereby making this approach valuable considering that 
exact computational schemes are very expensive. 



2. Mean-field theory 

The QLE description of Sec. IIII A II can be generalized to accommodate electron-electron repul- 
sion effects on the dot at the mean-field level. We refer to this extension as a MF QLE treatment, 
and note that it is not trivial: While a MF theory has been developed, suffering from some patholo- 
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FIG. 2: (a) Dot occupation as a function of time calculated using either the QLE method (full line) or 
the trace formula (dashed line), described in Sec. IIIIBI (b) Population of two selected states of the left 
reservoir, plotted as a function of time, e/ = —0.2 (top) and e/ = (bottom). In both cases we simulated 
the noninteracting Anderson model with £d = 0, \Il = —fJ,R = 0.2, Tl = Tr = 0.025 and = 0, Nlji=201 
and D = ±1. 



gies, for the study of dot occupation or charge current in the steady-state limit [33, |38|], here we 
present a MF scheme to describe the real-time dynamics of the full density matrix. By comparing 
MF results to exact numerical simulations, see Sec. IIV[ we conclude that a MF description can 
produce physical results up to ^, flL U ^ llR ^5 0{1). The effectiveness of the method also delicately 
depends on the dot level position, see for example Fig. [71 

The MF prescription, treating Coulombic repulsion, takes us back to Eq. ([2]). We now assume 
that the many-body interaction term can be factorized 



37, 



38] 



Un dA n d ^ ->■ U [(n djt ) n dA + n d)t (n dti )] . 



(25) 



This assumption reduces the Hamiltonian to an effectively noninteracting one, with a renormalized 
dot energy 



e d (t) = e d + U (n d (t)) 



(26) 



The spin-index has been dropped here, as we choose not to study magnetic effects. The formalism 
could be feasibly generalized to include magnetic fields, resulting in e^-f- ^ e d ^. In such situations 
the validity of MF equations is governed by another energy scale besides U/T and U/Afi, namely 
U/(e d + — e d) i). The dot occupation is determined in a self-consistent manner at every instant by 
modifying Eq. (|19p to contain the dot renormalized energy, 



(Mt)) = £ 



\v k \ 2 fk 



l + e 



-2Tt 



2e- it cos[(e dfc + U(n d (t)))t] 



(27) 



The solution provides the renormalized dot energy i d (t), which is then used to replace e d in Eqs. 
SB), flM}, ([A3]) and (dU), to provide the full DM at the MF level. 
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3. Probe model 



To implement elastic dephasing or inelastic effects with a probe, the set of equations (|12|) is 
augmented by an additional equation for c g . The EOM for a must be modified to include its 
coupling to the G bath, 

Cg ItgCg IVgCrf 

Cd = -ierfCrf -i}, vici - i/, v r Cr - i VgCg- (28) 

l r g 

It can be easily shown that under the QLE basic assumptions, as discussed in Sec. IIII A 11 the dot 
still satisfies Eq. (|17p with an additional noise term rp and with a re-defined total hybridization, 
T = ri+r^ + rc. The noise rp obeys a relation analogous to Eq. (fT5j) . and Eq. (|18p is generalized 
to 

Cd (t) = c d (t )e ( - 4£d - r)( *-* o) - i [ e H^-r)(*-r) [t/ l (t) + ^ (r) + ^G( T )] dT . (29) 

Jt 

We can now recognize that, in the presence of the probe, the expressions for the DM elements (1191) . 
(|20|) . (j'24|) . ()A3p and (|A4p stay formally intact. The technical adjustments are as follows: (i) We 
re-define the total hybridization, T = Tl + + ^g- (h) We augment summations that run over 
both L and R baths by k € G terms. For example, the summation in F% [Eq. (|A2|) ] should include 
such terms, (iii) We set the G bath distribution to satisfy the probe conditions, explained in Sec. 
d 



B. Fermionic trace formula (U = 0) 

We describe here an exact brute force calculation that can provide numerically all the elements 
of the density matrix in the noninteracting case. We begin without the presence of a probe, opting 
to include its effects later. This unitary method complements the QLE description, whose validity 
is governed by r^. Since the method is unitary, a recurrences behavior is ex pec ted to manifest at 
long enough time. The core of the method is the trace formula for fermions 13J 

Tr [e Ml e Ma ...e M *] = det [1 + e mi e m2 ...e m "] , (30) 

where m p is a single-particle operator corresponding to a quadratic operator M p = ^ j( m p)i,j c \ c j- 
c] (cj) are fermionic creation (annihilation) operators. Our objective is the dynamics of a quadratic 
operator A, either given by system or bath degrees of freedom, A = cjc&, j, k = l,r,d, 

(A(t)) = Tr [p(t )e iHt Ae- iHt ] = lim^o^Tr [p L PRP d e lHt e XA e- im ] . (31) 

13 



We introduce the A parameter, taken to vanish at the end of the calculation. The initial condition 
is taken to be factorized, p(to) = Pd® PL® PR, Pv = e~^^ H "~^" N ^ /Z v , Z v is the partition function, 
Pd describes the dot initial density matrix. These density operators follow an exponential form, 
e M , with M a quadratic operator. The application of the trace formula leads to 

(e^W } = det { [I L - f L ] ® [I R - f R ] [I d - f d ] + e M e Xa e - iht f L ®f R ®f d }. (32) 

Here, a and h are single-body matrices of the A and H operators, respectively. The matrices I v 
and Id are the identity matrices for the v = L,R space and for the dot. The functions fi and 
fn are the band electrons occupancy /^(e) = [e^ e ~^ v ' + 1]~ 1 . Here they are written in matrix 
form and in the energy representation, fd represents the dot initial occupation, again written in 
a matrix form. Since we are working with finite-size reservoirs, Eq. (|32p can be readily simulated 
numerically-exactly. 

The fermionic trace formula can be trivially generalized to include a probe. We add the G bath 
into the expectation value expression, A = cjcfc, 

{A(t)) = Tr [p(t y Hpt Ae- iHpt ] = lim A ^Tr [p LP RPGP d e iHpt e XA e- iHpt ] , (33) 

where as before p v = e ~^{ H "—^ N ^) jz v ^ Z v is the partition function, v = L,R, G. p d stands for the 
dot initial density matrix, and we trace over all DOF, the two reservoirs, the probe, and the dot. 

Timescales. Simulations with the trace formula are not restricted to a certain time scale. The 
method is unitary, providing (physical) recurrence behavior due to finite size effects. Since the 
time evolution scheme is not iterative, the accuracy of results does not deteriorate in time. 

C. Numerically exact path integral simulations, U ^ 

The time evolution of the closed and interacting Anderson model can be simulated by em- 
ploying a numerically-exact iterative influence- functional path integral (INFPI) approach [lj, Q]- 
This method relies on the fact that in out-of-equilibrium (and nonzero temperature) cases bath 
correlations have a finite range, allowing for their truncation beyond a memory time dictated 
by the voltage-bias and the temperature. Based on this finite-memory assumption, an iterative- 
deterministic time-evolution scheme has been developed, where convergence with respect to the 
memory length can, in principle, be reached. The principles of the INFPI approach have been 



detailed in Refs. 
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la], where it has been developed to investigate dissipation effects in the 
nonequilibrium spin-fermion model, and the population and current dynamics in correlated quan- 
tum dots. Recently, it has been used to examine the effects of a magnetic flux on the intrinsic 

14 



coherence dynamics in an Aharonov-Bohm quantum dot interferometer [39|. The INFPI method 
relies on the existence of a finite decorrelation time, thus it is suited for simulating the dynamics 
of an impurity coupled to a bath. Here we show that it can be used to retrieve the total DM in a 
system-bath setup. While in principle the method could encompass both interactions and probe, 
we focus exclusively on the first element. 

The method is based on the fermionic trace formula (|30p . incorporating many-body effects 
within a path integral expression. Our work starts with the time evolution expression (|3ip under 
the Hamiltonian ([2]). We factorize the time evolution operator, e lHt = (e lHSt ) Nt , NfSt = t, and 
adopt the Trotter decomposition e iH6t ^iHoSt/2 e iHiSt e iH st/2^ w h ere H = H + H x with 

u=L,R - ^ 



Hi = U[n d:t n d ^ - -(n djt + n d ^)] . (34) 



2 

H\ extracts many-body interactions on the dot, and it is eliminated 
variables s = ± via the Hubbard-Stratonovich (HS) transformation 



introducing auxiliary Ising 



2 



Here, k± = K'^fin", k' = sinh _1 [sin((5i[//2)] 1 ' /2 , k" = sin -1 [sin(<5£[//2)] 1//2 . The uniqueness of 
this transformation requires that U5t < tt. Incorporating the Trotter decomposition and the HS 
transformation into Eq. (|31|) . we find that the time evolution of A is dictated by 

(A(t)) = lim^{ j dsfds±..ds% t I(sf,st,...,s± t )}. (36) 
The integrand, referred to as as the "Influence Functional" (IF), is given by (q = 1, q + p = Nt) 

(37) 

where £+(s+) = ( e iH St/2 e H + (s+) e iH St/2\ and g_ = Eq (J36J) is exact in the St -> limit. 



Practically, it is evaluated by truncating the IF beyond a memory time r c = N s 5t, corresponding 
to the time beyond which bath correlations may be ignored N s is an integer. The following 



(non-unique) breakup has been suggested by [15j |. 



I( s fi s 2 = j ••• s jVt) — I( s f> s 2i - i • S W S )^( S 2 = ' s 3' s ]v s +l)"'^ s ( S iVt-Ar a +ls s N t -N s +2 J ■ » S Ar t )' 



where each element in the product, besides the first one, is given by a ratio between truncated IFs, 

^X 5 ^' S <j+1' ••> S q 
-^( S 9 J S g+1' •■) S q+i 



Tt \ ^( s ?' S 9+i'" -,S 9+Ar s -i) , QfV . 
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We now define a multi-time object, 



s ± s ± 



and evolve it iteratively by multiplication with the subsequent truncated IF, followed by summation 
over the time variables at the head, 

^■( s g+2> S g+3' •"' S q+N s ) = ^( s q+l' S q+2' — > S g+AT s -l)^( S g+l' s q+2' "•' S q+ivJ- (^1) 

The behavior at a particular time t q is reached by summation over the internal variables, 

{e XA(t q))= £ n(sf +2 _ Ns ,sf +3 _ Ns ,...,s±). (42) 

s ± s ± 
S q+2-N 3 > — > s Q 

This procedure is repeated for several (small) values of A, and the expectation value (A(t q )) is re- 
trieved by numerical differentiation in A. The truncated IF, Eq. (|37j) . is the core of this calculation. 
It is achieved numerically-exactly using the fermionic trace formula (|30|) . 

Timescale. Previous studies for dense reservoirs have confirmed that INFPI can provide accurate 
^ in both sh o rt tim e m d in th. ^i stea^e reg io„ QQ. However, t „e m e t hod is n „ t 
restricted to such dense-reservoirs situations, and it can describe the dynamics of small metallic 
grains since it handles all states explicitly. It should be still noted that the basic working assumption 
behind INFPI is the existence of a finite bath-induced decorrelation time. If the metal grains are 
very small, including few discrete states, this memory time r c does not exist or it becomes large, 
hindering convergence. Roughly, one could expect that a decorrelation time can be identified when 
a system-bath picture still holds, in the sense that a QLE description can be written i.e., Eq. (|17|) 
is valid. In such situations, INFPI simulations should converge and generally hold beyond t&. In 
practice, since these calculations are intensive, we have computed dynamics within a relatively 
short interval, Tt < 5, where the QSS description is still valid. 

IV. APPLICATIONS 

We now turn our attention to applications of the preceding methods. We first study the effects 
of Coulombic interactions on the reservoirs' DOF evolution. We later investigate the equilibration 
process in the system, through the implementation of probes. 
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FIG. 3: Population of reservoirs' levels in the noninteracting case with the dot energy positioned (a) within 
the bias window at = 0, (b) above the bias window at td = 0.3, and (c) below it at e<j = —0.3. Plotted arc 
the L (three top lines) and R (three bottom lines) occupations as a function of electron energy, at Tt = 
(dashed) Tt — 9.5 (full) and Tt — 19 (dotted). The framework used is a quantum Langcvin approach (Sec. 
IIII Aj) with (3 = 200 for the inverse temperature, Tl = Tr = 0.025 for bath-dot hybridization, r = Tl + Tr, 
fj,L = —fiR = 0.2 as a symmetric bias. The reservoirs are modeled by flat bands with a sharp cutoff at 
D = ±1, including N = 501 electronic states for each reservoir. 
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FIG. 4: Population of reservoirs' states from a mean-field QLE treatment. Plotted are the L (three top 
lines) and R (three bottom lines) occupations as a function of electron energy at Tt = 15 using the same 
set of parameters as in Fig. [31 (a) = 0, (b) = 0.3, and (c) td = —0.3. Full, dashed, and dashed-dottcd 
lines correspond to U = 0, 0.1, 0.3, respectively. The initial distributions for both the L and R reservoirs are 
presented by dotted lines. 



A. Anderson model with electron-electron interactions 



In this section we study the evolution of the finite-size Anderson model with or without in- 
teractions, based on the three methods described earlier in Sec. IIHI As mentioned above, these 



17 



techniques provide the dynamics of the total DM. While the fingerprints of many-body effects 
are disguised in the time evolution of conventional quantities, e.g., in the dot occupation and the 
charge current, they are well manifested in the reservoirs' population dynamics, allowing us to 
discern microscopic many-body scattering processes from single-particle events. 

The population p(e^) = (c^c^) of both reservoirs, in the noninteracting case, is displayed in Fig. 
[3]at different times. We note that results obtained using the QLE framework of Sec. IIIIAI perfectly 
agree with numerically-exact fermionic trace formula simulations. The three panels present results 
using different values for the dot energy, (a) When the dot energy is placed within the bias-window 
(e<i = 0) a resonance feature develops around the position of the dot level, with a dip (peak) showing 
in the L (R) bath. In contrast, if the dot energy is positioned either above the bias window (b) 
or below it (c), a dot-assisted tunneling feature develops, with population transfer taking place 
around available states that are the nearest in energy to e^. The dynamics shown in Fig. [3] is 
reversible, with a characteristic time ~ 2n/AE, AE = 2D/N is the mean spacing between 



361 ] . At this characteristic 



energy levels and N = N^r is the number of states in the L, R baths 
time the dot population begins to vary from its QSS value due to finite size effects. This behavior 
can be captured with trace formula simulations, but not within the QLE approach. 

In Fig. H|we display the dynamics under a mean-field QLE treatment with parameters corre- 



381 ] , we stress that 



sponding to Fig. While we are mindful of the technique's known pathologies 
this calculation provides an intuitive understanding of the role of interactions: Within MF, the ef- 
fect of finite U is to shift features in concert with the renormalized dot energy, edit) = e^ + U (rid(t)) 
[Eq. (I26p ]. Interestingly, panel (c) demonstrates a change in transport mechanism, from a dot- 
assisted tunneling at small U, to resonance transmission at large U, since the renormalized dot 
energy enters the bias window at a large enough interaction strength. Therefore, e-e interactions 
can enhance or suppress electronic transport, depending on the dot bare energy position. 

Mean- field results are compared to numerically-exact INFPI simulations in Fig. [5] for f7 = 0.1 
and U = 0.3, with the bare dot energy centered within the bias window. Data was produced by 
time evolution of all (clc^), k = l,r, expectation values up to Tt ~ 4. In agreement with MF 
QLE results, the basic effect of e-e interactions observed here is a shift in the resonance position. 
Overall, we conclude that MF simulations can reproduce the dynamics for this set of parameters, 
up to U/T < 2. Qualitative features are correct through U/T ~ 6. 

Convergence of INFPI is verified with respect to the time step adopted, 6t, and the memory 
time accounted for, r c . Representative convergence curves for p(ei = e^) are depicted in Fig. 
[6l While the time-step used does not affect our results, we note that the data is not yet fully 
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converged with respect to r c . This slow convergence could be attributed to the long decorrelation 
time experienced by an electron residing on any particular bath level, since its decorrelation process 
should take place by following a two-step procedure: the electron should first leave the particular 
bath state and populate the quantum dot. From the dot, it may subsequently transfer to any other 
bath state. One should also note that we display here a convergence curve for a particular level. 
It is more accurate, but computationally demanding, to look at the overall evolution of p(e) (for 
all e) with r c . This is because the resonant feature may change its magnitude with r c , as we see 
here, as well as its absolute position. Fig. [6] only analyzes the effect of r c on the peak magnitude. 
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FIG. 5: Simulations of reservoirs' levels population in the resonance case, = 0, at Tt = 4, using INFPI 
(dark curves) and MF QLE (light curves). Full, dashed and dashcd-dottcd lines were obtained using U = 
0,0.1,0.3, respectively. Top lines correspond to the L bath distribution, the bottom lines to the R bath. 
The initial distributions of both L and R reservoirs are presented by dotted lines. Parameters are the same 
as in Fig. [3J with Nl r = 101 bath states. INFPI numerical parameters are St = 1 and N s = 7. 
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FIG. 6: Convergence trend for the INFPI data of Fig. [5] We plot p(ei = e^) as a function of r c at Ft = 1.2 
for U = 0.1 and U = 0.3. Results are shown using different time steps, St = 0.8 (□), St = 1.0 (V), St = 1.2 
(0), St = 1.5 (x) and St = 1.6 (o). Other parameters are the same as in Fig. El 
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FIG. 7: (a) INFPI data for the reservoirs population at different times for = 0.25 and U = 0.1. The 
five top (bottom) lines depict the population of the L (R) baths at times Ft = 0.05, 1, 2, 3, and 4. (b) 
Zooming over the L bath population, (c) Zooming over the R bath population. The arrows indicate on the 
direction of time evolution. In (a)-(c) the full lines represent MF QLE data at the latest time, Ft = 4. (d) 
Convergence behavior of p(e r = e^) at Ft ~ 1.2 using different time steps, St = 0.8 (□), St = 1 (v) and 
St = 1.5 (o). Other parameters are the same as in Fig. El with Nl^r = 101 bath states. 



In Fig. [7] we study the dynamics of the reservoirs' levels population when the renormalized 
dot energy sits above the bias window, > /i^. The bare energy is taken at = 0.25 and the 
interaction strength is U = 0.1. By separately calculating the time evolution of the dot occupation, 
to produce (n^a) ~ 0.12 in the QSS limit (valid for Tt > 2), we estimate the renormalized dot 
energy to be about ~ 0.26. Overall, occupations change very slightly in time, since the dot is off- 
resonant thus the transport follows a dot-assisted tunneling mechanism. The reservoirs' dynamics 
still clearly manifests many-body effects that are not included in a MF (effective single-body) 
description. Essentially, we find that electrons populate high energy levels in the L and R baths 
up to ek ~ 0.45. This high-energy population cannot be explained by the dot-level shift or the 
dot finite broadening T, as this could account for population up to e& ~ 0.35 only, see MF data 
(full line) in panel (c). It should be noted that the population of levels that are initially empty, 
p(e > develops identically at the L and R reservoirs. In other words, the high energy tails in 
panel (c) represent the occupation of states both in the L and R baths. Supporting convergence 
behavior is included in panel (d). We have also performed simulations when taking a stronger 
interaction, U = 0.3 and = 0.15, to yield e ~ 0.21. In this case the population shows a high- 
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energy tail that is further enhanced with respect to MF data, representing significant deviations 
from a single particle description. However, as we did not manage to fully converge these results, 
they are not included here. 

The breakdown of the single-particle picture is difficult to discern in cumulative quantities such 
as the charge current and the dot occupation, the latter is presented in Fig. [SJ where we follow 
the time evolution of this quantity using four different techniques: MF QLE equations, first-order 



perturbation theory method 



341 ] . Monte-Carlo simulations 
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44 1 and INFPI [131 . The comparison 

n 



shows a very good agreement between the latter two exact methods up to Ft ~ 1 [151 ]. MF QLE 
theory and first-order perturbation expansion both predict the correct behavior at U/T < 2. At 
strong interactions (U/T ~ 6), perturbation theory fails, while MF QLE equations still provides 
correct qualitative behavior for (rid t a)- Note that we implement a sharp cutoff at D = ±1 in all 
methods. Since we do not account for the (small) energy shift in Eq. (|17p . MF QLE data suffer 
from a small shift in values, see also Fig. [2] 

Naively considering the dot occupation only, one may conclude that at U/T < 6 many-body 
effects are contained in a mean-field, effective single-body, description. However, traces of energy 
resolved reservoirs' dynamics as we show in Fig. [7]expose the existence of interaction effects beyond 
mean-field, resulting in levels population beyond the resonance width. The dot occupation thus 
withholds mechanisms involved in the transport process, while detailed reservoirs' level population 
can illuminate extant many-body effects and their energy resolution. 



U=0.3 




FIG. 8: Dot occupation as a function of time, generated from four different methods: INFPI (dotted), 



Monte-Carlo simulations (□) |4lU44j). MF QLE equations (dashed) and perturbation theory treatment (o) 



34 1 . Results are shown for three setups, bottom to top: = 0.3 with U = 0, = 0.25 with U = 0.1, and 
£d = 0.15 with U = 0.3. Other parameters are the same as in Fig. [3] with 8t = 1, N s = 7 and Nl,r = 101 
bath states. 
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We now comment on the simulation time of a convergence analysis as presented in Fig. [6l 
covering three different time steps and N s = 2, 7. Convergence should in principle be verified for 
all bath states (Nl+Nr = 101 x 2) and at all times. In practice, we have tested it for representative 
states only. The MATLAB implementation of the computational algorithm took advantage of the 
MATLAB built-in multi-threaded parallel features and utilized 100% of all available CPU cores on 
a node. When executed on one cluster node with two quad-core 2.2GHz AMD Opteron cpus and 
16GB memory, convergence analysis for each expectation value took about 7x24 hours and 250MB 
of memory. Computations performed on the GPC supercomputer at the SciNet HPC Consortium 



45J were three times faster. Computational time scales linearly with the simulated time t. For a 
fixed N s value, the computational effort does not depend on the system temperature and the value 
of U employed. 

B. Quantum equilibration and thermalization 

The techniques developed in Sec. IIIII provide the time evolution of the total DM of a large 
system, allowing us to address next the problem of equilibration and thermalization in quantum 
mechanics. The basic question of interest here is how do quantum systems equilibrate from a 
nonequilibrium initial preparation, if at all. Furthermore, it is of importance to understand under 
what conditions a system may approach one of the canonical ensembles of statistical mechanics. 

Before addressing the equilibration problem in detail we present in Fig. [9] a more standard 
quantity, the steady-state dot occupation. This will serve us for motivating the study of the total 
DM for resolving transport mechanisms. The dot occupation is displayed here as a function of 
Tq, the dot-probe hybridization strength, and we show results using either a voltage probe or a 
dephasing probe, for different values of the dot energy position. We find that the dot occupation 
is insensitive to the probe condition, an observation that can be proved analytically by studying 
the long time behavior of Eq. (|19j) under either probes, Eqs. fjlQf) or (llip . It is interesting to note 
the crossover behavior of the dot occupation when its energy is placed above the bias window: 
When Tg < ^l,r occupation grows linearly with Tq. However, it decays as T^ 1 at large values, 
when effective dephasing and inelastic effects are strong. This behavior is similar to the thermally 
assisted tunneling behavior observed when using more detailed modeling [46J. It can similarly be 



shown that the steady-state current in the system is identical irrespective of the probe condition. 
Note that we restrict ourselves to the quasi steady-state region since the G bath does not serve as 
a proper probe before quasi steady-state sets in. 
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The underlying transport mechanisms are therefore obscured in cumulative quantities (current, 
occupation) in this probe model, as well as in more explicit electron-phonon modeling j^]. Details 
about the involved mechanisms can be resolved by studying, e.g. the current noise, inelastic electron 
tunneling spectra 



471 ]. and the evolution of the reservoirs population 



481 ] . as we show below in the 



context of quantum equilibration and thermalization. 



0.8 
0.6 
0.4 
0.2 




\e d =-0.3 



e=0 --O-- 



a 







0.5 



1.5 



FIG. 9: Steady-state dot occupation under a voltage probe condition (dashed) and a dephasing probe (dot) 
for different values of the dot energy, = —0.3, and 0.3, top to bottom. Data was obtained by simulating 
Eq. (TJU) in the long time limit using fi L = -/j, R = 0.2, (3 = 200 and T L , R = 0.025. 



The equilibration 



problem in quantum mechanics could be considered within different setups: a 



closed system 49to1[] . a system-bath scenario 52J, [53J], or taking peer quantum systems 25f] . Here, 
we consider the Anderson model with a probe, excluding e-e interactions, and simulate the following 
setup: At time t = to we put into contact through a quantum dot two reservoirs each separately 
prepared in grand canonical states at chemical potentials fi u and temperature /3 _1 . Electrons on the 
dot are susceptible to either elastic-decoherring processes or inelastic effects, mimicked by coupling 
to the relevant probe. For a schematic representation, see Fig. QJb). Given this scenario, we 
investigate whether the two reservoirs can equilibrate or even thermalize in time, and furthermore, 
the nature of the equilibrium state. As we show below, when only elastic dephasing effects are 
mimicked, the system approaches a non-canonical equilibrium state. When inelastic processes are 
emulated, the two reservoirs relax towards a common canonical state. It should be noted that 
these results can be obtained for a finite and closed system, under a unitary evolution 23]. 



We identif y th ermal equilibration in our peer quantum system setup, by adjusting the conditions 
of Refs. 



52, 



54j , demanding that: (i) The system should equilibrate, i.e., evolve towards some 



particular state, and stay close to it for almost all time. Furthermore, the equilibrium state should 
be (ii) independent of the dot properties-energetics and initial state, (iii) insensitive to the precise 
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FIG. 10: Dynamics of the reservoirs population in the Anderson probe model, (a) Population of the L 
(full) and R (dashed) reservoirs, Tq — 0. We show results at t = 0, 180, 340 as the resonance pattern 
develops in time, (b)-(c) Equilibration with a dephasing probe, Tg = 0.4, approaching a non-canonical 
equilibrium state. The L (full) and R (dashed) population are shown at t = 0, St, 9St, St = 600. Panel (c) 
demonstrates a slow-down in dynamics in approaching the equilibrium state, (d)-(f) Approaching thermal 
equilibrium with a voltage probe. In panel (d) Tq ~ 0.4 with the L (full) and R dashed line population 
shown at t = 0, 700, 2800, 6300. Panel (e) displays the population as a function of time at a certain energy, 
e = 0.1. Panel (f) presents information as in (d), with = 5. In all panels /3 = 200, Tl = = 0.025, 
e c i = 0, /.iL = —fJ-R = 0.2, D = 1, Nl.r = 101 and Ng = 2001. The arrows mark the direction of propagation 
in time. 

initial state of each reservoir, (iv) close to diagonal in the energy basis of its eigen-Hamiltonian, 
and (v) a canonical state. 

We use the trace formula, an exact unitary method, and follow the reservoirs' mutual equilibra- 
tion process. We evolve the system using either a dephasing probe or a voltage probe, see Fig. [TOl up 
to the time where recurrence features start to manifest, found here to scale as r rec oc ^2 i=L fjgiV,. 
As a reference, panel (a) displays results for the model without a probe, showing the development 
of a resonance feature around the dot energy position at = 0. A clear evolution towards an 
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equilibrium state is demonstrated when a probe is presented. With a dephasing probe, (b)-(c), 
the population of the two reservoirs relax to a two-step function with p(fj,R < e < \xl) ~ 0.5. 
Because electrons from the L grain loose their phase memory on the dot, half populate the R 
side, on average, in the long time limit. This equilibrium state is sensitive to the precise details 
of the initial electron distribution, as energy redistribution is not allowed. We build a large G to 
delay recurrence behavior, but note that results at earlier times do not depend on the size of G, 
reinforcing the observation that G acts as an agent in driving the L-R mutual equilibration. When 
inelastic effects are mimicked with a voltage probe, and Tq is large enough, panel (f), the system 
approaches a Gibbs-like thermal state — a step function at zero temperature. Results are shown up 
to the time r rec at which recurrence features develop, which emerges here before full thermalization 
takes place. In order to achieve full thermalization one should further increase the size of the G 
bath, so as to delay recurrences. Alternatively, a dissipative mechanism could be introduced into 
G, e.g. by building a hierarchy of its interactions with the L-R system. Using a smaller value 
for Tq, a non-canonical equilibrium distribution develops (d)-(e), reflecting the contribution of 
coherent and (effectively) incoherent electrons in the dynamics. It is also interesting to compare 
panels (c) and (e), displaying the equilibration progress for a dephasing probe and a voltage probe, 
respectively, while maintaining the value of Tq. We find that the characteristic timescale to reach 
equilibrium is very similar in both cases. Thus, while the probe type dictates the structure of the 
equilibrium state, it does not affect the equilibration timescale. 

Fig. QT]shows that while under coherent evolution the resonance peak emerges around the energy 
€d, in the presence of a voltage probe with (large enough) Tq, the buildup of the equilibrium state 
systematically occurs around the equilibrium Fermi energy. This holds even when the dot is placed 
outside the bias window (not shown). Analogous trends take place when allowing for dephasing 
only. 

A thermal equilibrium state should be diagonal in the energy eigenbasis of its Hamiltonian 521 ]. 



In Fig. [12] we display the density matrix = (etc*;/), k,k' = l,r, excluding diagonal elements 
Pk,k, without a voltage probe (a)-(b), and with a one (c)-(d), using the QLE technique. This 
quantity is expected to oscillate in the long time limit since the Hamiltonian is not diagonal in the 
(local) / and r bases. We still show the results in the local reservoirs' basis, so as to manifest local 
z^-bath properties. There are three significant differences in the behavior of off-diagonal elements, 
with and without the probe: (i) The absolute value of the coherences, at a given time, is smaller 
when Tq ^ 0. (ii) The DM approaches a diagonal form (strict diagonal values are not shown), (in) 
When Tq = 0, oscillations occur around the dot energy position e^. With the probe, contributions 



25 




-0.4 -0.2 0.2 -0-4 -0.2 0.2 



FIG. 11: (a) Occupation of L bath at time t = (o) and t = 1500 for T G = and e d = (dotted), T G = 
and ed = 0.1 (dashcd-dottcd), T G = 0.4 and e d = (full), T G = 0.4 and e d = 0.1 (dashed). The latter two 
lines assume the voltage probe condition, (b) Same for the R side occupations. Other parameters are the 
same as in Fig. [TO] 

are scattered, yet they appear more prominently around the equilibrium Fermi energy, e = 0. These 
three features should become more pronounced at longer times, which can be simulated using the 
trace formula approach. 




FIG. 12: Absolute values of the density matrix elements Pk,k', k,k' = l,r, at t = 600, excluding diagonal 
elements. In panel (a) and (b) Tq = 0. In panels (c) and (d) we use a voltage probe with Tg = 0.4. Panels 
(a) and (c) display the total density matrix, and the axes are the energy indices 1,...,402. The L reservoir 
includes the first 201 states. The rest are R bath states. Panels (b) and (d) zoom on the density 
matrix, the bottom-leftmost part of the total DM. Other parameters are ed = 0.1, = 0.025, j3 = 200, 
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V. CONCLUSION 



We have extended analytical and numerical methods, developed for simulating the dynamics of 
impurities, i.e., subsystems attached to large reservoirs, to reveal the dynamics of the total system. 
As an example, we have focused on the Anderson model, a quantum dot coupled to two metal 
grains, and obtained the evolution of the total density matrix, focusing on the reservoirs' evolution 
from an initial nonequilibrium state. We have studied the noninteracting model, as well as a 
model with interactions and a probe model, emulating elastic dephasing and dissipation effects. 
The three methods presented are the analytic quantum Langevin equation approach, a simulation 
based on a trace formula, and an exact numerical path integral scheme that can accommodate 
e-e repulsion effects. Notably, the extension of the QLE treatment to provide the total DM is of 
general importance as it can be used in multitude of other systems, as long as one can identify a 
"subsystem" within the total system. 

Making use of the methods developed, we have investigated the total system dynamics in the 
presence of distinct effects: (i) e-e interactions on the impurity, and (ii) dephasing and inelastic 
scattering effects. Addressing the prior, our calculations allow us to energy resolve the effect of 
e-e interactions on electron transfer in the Anderson dot model. In the resonant regime we found 
that the dynamics observed for noninteracting electrons is largely preserved up to U/Y < 2, and 
the main effect of interactions on the reservoirs' occupation is apparently a simple shift in the 
position of features affected by the renormalization of the dot energy. Away from resonance, in the 
tunneling domain, the presence of weak interactions already manifested itself in scattering electrons 
to high energy levels, an effect that is not captured within a mean- field treatment. In the case 
of the later effect, we found numerically that the presence of dephasing and inelastic effects on a 
weak link only can lead to global system equilibration and even thermalization. It is important to 
note that no restrictions were enforced on the metals' band structure and the dot energy. This is 
significant in light of many other studies in which equilibration requires the "nondegenerate energy 



gap" condition to be satisfied [25 



49, 



Future directions include the study of finite temperature and electron-electron interaction effects 
in the equilibration process and the behavior given a quantum dot chain between the two metal 
grains. In a linear _chain of impurities we expect that the coherent-diffusive crossover in the charge 
current behavior 



55J would similarly manifest itself in the energy reorganization process of the 
reservoirs. The methods developed here could also be adopted for the study of bosonic systems, 
e.g., to describe the dynamics of bosonic degrees of freedom interacting with harmonic baths. 
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Appendix A: Density matrix elements in the quantum Langevin approach 



We provide here the explicit expressions for the density matrix elements pik, k € L,R. The 
population is calculated by evaluating the correlation functions in Eq. (|24p . using Eqs. (|13p and 
(USD, to yield 



p(ei) = {cJ(t )ci(t )) + F 2 + F 3 



(Al) 



The first term accommodates the initial condition. The second (F2) and the third (F3) terms are 
given by 



F 2 = -vffi 



2T 



r2 + 



-t-2vff l 



di 



:4 + r 2 ) 



+ 



vffie 



-17 



dl 



2 (4 - T 2 ) cos {e dl t) + 4e^r sin (e dl t) 



F 3 = v? 



v k 2 f k J 4 sin 2 (^) 1 



r 2 + 



+ 



+ 



I _ e -K r + i€ dl) _|_ e -K r + ie dk) _ e ~ ite 



r 2 + 



di 



e -2Ft + 1 _ e t(ie dl -r) _ e -t(ie dl +T) 



+ c.c. 



(A2) 



(e d i - iT) eik 

Inter and intra-reservoir coherences, e.g. pi ik (t), k € L,R can be similarly calculated. Here, one 
should distinguish between the cases €/ = e k and e\ ^ e k - In the latter case we find that 

A,fc(t) = (4(*) c fe(*)) = A X +A 2 + A 3 (A3) 



where 



At 



A-, 



VkVlfl 

r + ie d i 

VlVkfk 

r - ie dk 



-t(T+ie d i) _ pieikt 



r + ie dk 

-t(F-ie dk ) _ ie tk t 



T - ie d i 



ezfc 



2S 



and 



-43 = VlVk 



V V k'fk' 

1 e ite ifc' e iteik e ite k'k e ite k'k + e -*(r+iedi) _ g«te;fc _ e -t(r+«e dfe /) 

+ + 



^k'l^k'k ^lk' e k'k Clk'tk'k ^k'l^k'k 



g* Ie ifc' _1_ e -*(r-iedfe) _ giteifc _ e -t(T-ie dk i) e ~ 2rt -L e ite lk _ g-*(r+ie di ) _ g-t(r— ie dii ) 

efc/fc(ir + e,a) («r + e di )(-«r + e d k) 

In the resonant limit, ei = and k ^ L, a simpler result is obtained, 

c](t)ck(t))=Al+A r 2 + Al 



(A4) 



with 



A r 2 



VkVlfl 
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